Friction in Constant Area Duct

Table of Contents

1. Introduction
2. Problem Workthrough
3. Appendix: Derivation

Introduction

In this notebook we will study how friction changes the Mach number in a constant area duct \ As we will see the coefficient of friction plays the same role as the change in area of a converging frictionless duct, meaning if the flow is originally subsonic the Mach number will increase and if it is originally supersonic the Mach number will decrease

# import relevant libraries
import plotly.graph_objects as go
import numpy as np
from scipy.optimize import root_scalar
from IPython.display import Markdown

Problem Workthrough

# Define constants and known values
M_1 = 0.5
x_1 = 2
c_f = 0.0025
D = 0.1
g = 1.4

If the pipe continues onto infinity at some point the flow will choke (reach M=1), this is useful because there is an isentropic flow relation between the Mach number and the distance to the choked position:

$$ \large \require{color}\frac{4c_{f}}{D}L = \frac{1 - {\color[rgb]{0.041893,0.355669,0.727621}M}_{1}^{2}}{\gamma {\color[rgb]{0.041893,0.355669,0.727621}M}_1^2} + \frac{\gamma +1}{2\gamma} \ ln\left(\frac{\frac{\gamma +1}{2}{\color[rgb]{0.041893,0.355669,0.727621}M}_1^2}{1+\frac{\gamma -1}{2}{\color[rgb]{0.041893,0.355669,0.727621}M}_1^2}\right) $$

This finds the maximum required distance for the flow to choke, then since we know the distance between $\require{color} {\color[rgb]{0.041893,0.355669,0.727621}M}_{1}$ and $\require{color} {\color[rgb]{0.041893,0.355669,0.727621}M}_2$ (where we want to know the Mach number) we can also find the distance it takes for the flow to choke from $\require{color} {\color[rgb]{0.041893,0.355669,0.727621}M}_2$ \ Then using the distance bwetween $\require{color} {\color[rgb]{0.041893,0.355669,0.727621}M}_2$ and $\require{color} {\color[rgb]{0.041893,0.355669,0.727621}M}=1$ we can solve the above equation for Mach number to find $\require{color} {\color[rgb]{0.041893,0.355669,0.727621}M}_2$ \ However, the equation above does not have an analytical solution so a process such as curve fitting or root finding must be used, in this notebook I use root finding using scipy

# define a function that calculates how far it takes for the flow to choke given the initial mach number 
def max_L(M):
    L = (D/(4*c_f))*((1-M**2)/(1.4*M**2) + (g+1)/(2*g) * np.log(((M**2)*(g+1)/2)/(1+(M**2)*(g-1)/2)))
    return L

# define the length to where the flow chokes and a position we want to find the Mach number 
L_choke = max_L(M_1)
x_2 = x_1 + L_choke/2

x=np.linspace(0, x_1+L_choke+1, 100)

# plot a duct and show where each Mach value is 
fig = go.Figure()
fig.add_scatter(x=x, y=np.e**(-x) + 1, fill ='tonexty', fillcolor='rgba(133, 163, 201, 0.5)', line_color='blue')
fig.add_scatter(x=x, y=-np.e**(-x) - 1, fill='tozeroy', fillcolor='rgba(133, 163, 201, 0.5)', line_color='blue')
fig.add_vline(x=x_1, line_dash='dash', annotation_text=r"$\require{color} {\color[rgb]{0.041893,0.355669,0.727621}M}_{1}$")
fig.add_vline(x=x_1+L_choke, line_dash='dash', annotation_text=r"$\require{color} {\color[rgb]{0.041893,0.355669,0.727621}M}=1$")
fig.add_vline(x=x_2, line_dash='dash', annotation_text=r"$\require{color} {\color[rgb]{0.041893,0.355669,0.727621}M}_{2}$")
fig.update_layout(yaxis=dict(scaleanchor="x", scaleratio=1, showticklabels=False), xaxis=dict(showticklabels=False),
                  showlegend=False, plot_bgcolor='rgba(0, 0, 0, 0)')
fig.show()

# define the length from M_2 to where the flow chokes
L_2 = L_choke + x_1 - x_2

# define function to find Mach roots 
def mach_function(M, L):
    return (D/(4*c_f))*((1-M**2)/(1.4*M**2) + (g+1)/(2*g) * np.log(((M**2)*(g+1)/2)/(1+(M**2)*(g-1)/2))) - L

# using scipy find the root where L=L_2
M_2_root = root_scalar(mach_function, args=(L_2), x0=M_1, x1=1)
M_2 = M_2_root.root

Markdown(fr'The Mach number at a distance ${round(x_2 - x_1, 3)}m$ from $M_1$ the Mach value is ${round(M_2, 3)}$')

The Mach number at a distance $5.345m$ from $M_1$ the Mach value is $0.589$

Appendix: Derivation

fig = go.Figure()
fig.add_shape(type='circle', xref='x', yref='y', x0=0, y0=0, x1=1, y1=3)
fig.add_shape(type='circle', xref='x', yref='y', x0=3, y0=0, x1=4, y1=3, line_dash='dash')
fig.add_shape(type='line', x0=0.5, y0=3, x1=3.5, y1=3)
fig.add_shape(type='line', x0=0.5, y0=0, x1=3.5, y1=0)
fig.add_annotation(x=0.5, ax=-0.5, y=1.5, ay=1.5, xref='x', axref='x', yref='y', ayref='y', arrowhead=5, 
                   text=r"$p+\rho v^2$")
fig.add_annotation(x=3.5, ax=5.5, y=1.5, ay=1.5, xref='x', axref='x', yref='y', ayref='y', arrowhead=5, 
                   text=r"$p + \frac{dp}{dx}dx+\rho v^2 +\frac{d\rho v^2}{dx}dx$")
fig.add_annotation(x=1.5, ax=2.5, y=3.2, ay=3.2, xref='x', axref='x', yref='y', ayref='y', arrowhead=5, text=r"$\tau_w$")
fig.add_annotation(x=1.5, ax=2.5, y=-0.2, ay=-0.2, xref='x', axref='x', yref='y', ayref='y', arrowhead=5)
fig.update_layout(yaxis=dict(scaleanchor="x", scaleratio=1, showticklabels=False), xaxis=dict(showticklabels=False),
                  showlegend=False, plot_bgcolor='rgba(0, 0, 0, 0)')
fig.show()

The shear force is constant across the surface area of the pipe spanned by a small section of fluid $dx$ this area is $2\pi R dx = \pi Ddx$ where $D$ is the diameter of the pipe\ The pressure and velocity terms only act on the inlet and exit areas which is given by $\pi R^2 = \pi \frac{D^2}{4}$\ Firstly, we must sum the forces:

$$ \large \require{color}F=-{\color[rgb]{0.990308,0.800015,0.121270}\tau}_w\pi Dd x + {\color[rgb]{0.315209,0.728565,0.037706}p} +{\color[rgb]{0.814433,0.253157,0.091125}\rho} {\color[rgb]{0.059472,0.501943,0.998465}v}^2 -\left({\color[rgb]{0.315209,0.728565,0.037706}p} + \frac{d{\color[rgb]{0.315209,0.728565,0.037706}p}}{dx}dx + {\color[rgb]{0.814433,0.253157,0.091125}\rho} {\color[rgb]{0.059472,0.501943,0.998465}v}^2 + \frac{d {\color[rgb]{0.814433,0.253157,0.091125}\rho} {\color[rgb]{0.059472,0.501943,0.998465}v}^2}{dx}dx\right)\pi \frac{D^2}{4} =-\pi \frac{D^2}{4}dx\left(\frac{d{\color[rgb]{0.315209,0.728565,0.037706}p}}{dx} + \frac{{\color[rgb]{0.814433,0.253157,0.091125}\rho} {\color[rgb]{0.059472,0.501943,0.998465}v}^2}{dx}\right) - {\color[rgb]{0.990308,0.800015,0.121270}\tau}_w\pi Ddx $$$$ \large \require{color}-{\color[rgb]{0.990308,0.800015,0.121270}\tau}_w\pi D = \pi \frac{D^2}{4}\left(\frac{d{\color[rgb]{0.315209,0.728565,0.037706}p}}{dx} + \frac{d({\color[rgb]{0.814433,0.253157,0.091125}\rho} {\color[rgb]{0.059472,0.501943,0.998465}v}^2)}{dx}\right) \longrightarrow -\frac{4{\color[rgb]{0.990308,0.800015,0.121270}\tau}_w}{D} = \frac{d{\color[rgb]{0.315209,0.728565,0.037706}p}}{dx} + \frac{d{\color[rgb]{0.814433,0.253157,0.091125}\rho} {\color[rgb]{0.059472,0.501943,0.998465}v}^2}{dx} $$

If we assume the flow conditions: $\require{color}{\color[rgb]{0.814433,0.253157,0.091125}\rho} {\color[rgb]{0.059472,0.501943,0.998465}v} = const, {\color[rgb]{0.064095,0.501831,0.501977}h} + \frac{1}{2}{\color[rgb]{0.059472,0.501943,0.998465}v}^2 = const$ then $\require{color}\frac{d{\color[rgb]{0.814433,0.253157,0.091125}\rho} {\color[rgb]{0.059472,0.501943,0.998465}v}}{dx}=0$ and $\require{color}\frac{d{\color[rgb]{0.064095,0.501831,0.501977}h}}{dx} + {\color[rgb]{0.059472,0.501943,0.998465}v}\frac{d{\color[rgb]{0.059472,0.501943,0.998465}v}}{dx}=0$ \ Using the ideal gas law $\require{color}{\color[rgb]{0.315209,0.728565,0.037706}p}={\color[rgb]{0.814433,0.253157,0.091125}\rho} {\color[rgb]{0.986252,0.007236,0.027423}R}{\color[rgb]{0.121820,0.954406,0.966585}T} $, the equation for enthalpy $\require{color}{\color[rgb]{0.064095,0.501831,0.501977}h}={\color[rgb]{0.501958,0.501942,0.014744}c_p}{\color[rgb]{0.121820,0.954406,0.966585}T}$, that $\require{color}{\color[rgb]{0.986252,0.007236,0.027423}R}={\color[rgb]{0.501958,0.501942,0.014744}c_p}-c_v$, and $\require{color}\gamma = \frac{{\color[rgb]{0.501958,0.501942,0.014744}c_p}}{c_v}$ we get that $\require{color}{\color[rgb]{0.315209,0.728565,0.037706}p}= \frac{\gamma -1}{\gamma} {\color[rgb]{0.814433,0.253157,0.091125}\rho} {\color[rgb]{0.064095,0.501831,0.501977}h}$\ Additionally since we know that $\require{color}{\color[rgb]{0.989013,0.435749,0.811750}a}=\sqrt{\gamma {\color[rgb]{0.986252,0.007236,0.027423}R}{\color[rgb]{0.121820,0.954406,0.966585}T}}$ and $\require{color}{\color[rgb]{0.059472,0.501943,0.998465}v}={\color[rgb]{0.041893,0.355669,0.727621}M}{\color[rgb]{0.989013,0.435749,0.811750}a}$ we can solve for $\require{color}{\color[rgb]{0.814433,0.253157,0.091125}\rho} {\color[rgb]{0.059472,0.501943,0.998465}v}^2$:

$$ \large \require{color}{\color[rgb]{0.814433,0.253157,0.091125}\rho} {\color[rgb]{0.059472,0.501943,0.998465}v}^2={\color[rgb]{0.814433,0.253157,0.091125}\rho} {\color[rgb]{0.041893,0.355669,0.727621}M}^2{\color[rgb]{0.989013,0.435749,0.811750}a}^2={\color[rgb]{0.814433,0.253157,0.091125}\rho} \gamma {\color[rgb]{0.986252,0.007236,0.027423}R}{\color[rgb]{0.121820,0.954406,0.966585}T}{\color[rgb]{0.041893,0.355669,0.727621}M}^2 = {\color[rgb]{0.315209,0.728565,0.037706}p}\gamma {\color[rgb]{0.041893,0.355669,0.727621}M}^2 $$

To expand $\require{color}\frac{d{\color[rgb]{0.315209,0.728565,0.037706}p}}{dx} + \frac{d{\color[rgb]{0.814433,0.253157,0.091125}\rho} {\color[rgb]{0.059472,0.501943,0.998465}v}^2}{dx}$ we expand the derivative of $\require{color}{\color[rgb]{0.814433,0.253157,0.091125}\rho} {\color[rgb]{0.059472,0.501943,0.998465}v}^2$ into terms of ${\color[rgb]{0.059472,0.501943,0.998465}v}$ and ${\color[rgb]{0.814433,0.253157,0.091125}\rho} {\color[rgb]{0.059472,0.501943,0.998465}v}$ and we use the expression for $p$ to expand the first term:

$$ \large \require{color}\frac{d{\color[rgb]{0.315209,0.728565,0.037706}p}}{dx} + {\color[rgb]{0.814433,0.253157,0.091125}\rho} {\color[rgb]{0.059472,0.501943,0.998465}v}\frac{d{\color[rgb]{0.059472,0.501943,0.998465}v}}{dx} + {\color[rgb]{0.059472,0.501943,0.998465}v}\frac{d {\color[rgb]{0.814433,0.253157,0.091125}\rho} {\color[rgb]{0.059472,0.501943,0.998465}v}}{dx} = \frac{\gamma -1}{\gamma}\left({\color[rgb]{0.814433,0.253157,0.091125}\rho} \frac{d{\color[rgb]{0.064095,0.501831,0.501977}h}}{dx} + {\color[rgb]{0.064095,0.501831,0.501977}h}\frac{d{\color[rgb]{0.814433,0.253157,0.091125}\rho}}{dx}\right) + {\color[rgb]{0.814433,0.253157,0.091125}\rho} {\color[rgb]{0.059472,0.501943,0.998465}v} \frac{d{\color[rgb]{0.059472,0.501943,0.998465}v}}{dx} \frac{d{\color[rgb]{0.315209,0.728565,0.037706}p}}{dx} + {\color[rgb]{0.814433,0.253157,0.091125}\rho} {\color[rgb]{0.059472,0.501943,0.998465}v}\frac{d{\color[rgb]{0.059472,0.501943,0.998465}v}}{dx} + {\color[rgb]{0.059472,0.501943,0.998465}v}\frac{d{\color[rgb]{0.814433,0.253157,0.091125}\rho} {\color[rgb]{0.059472,0.501943,0.998465}v}}{dx} = \frac{\gamma -1}{\gamma}\left({\color[rgb]{0.814433,0.253157,0.091125}\rho} \frac{d{\color[rgb]{0.064095,0.501831,0.501977}h}}{dx} + {\color[rgb]{0.064095,0.501831,0.501977}h}\frac{d{\color[rgb]{0.814433,0.253157,0.091125}\rho}}{dx}\right) + {\color[rgb]{0.814433,0.253157,0.091125}\rho} {\color[rgb]{0.059472,0.501943,0.998465}v} \frac{d{\color[rgb]{0.059472,0.501943,0.998465}v}}{dx} $$

Using that $\require{color}\frac{d{\color[rgb]{0.064095,0.501831,0.501977}h}}{dx} + {\color[rgb]{0.059472,0.501943,0.998465}v}\frac{d{\color[rgb]{0.059472,0.501943,0.998465}v}}{dx}=0 \longrightarrow \frac{d{\color[rgb]{0.064095,0.501831,0.501977}h}}{dx} = -{\color[rgb]{0.059472,0.501943,0.998465}v}\frac{d{\color[rgb]{0.059472,0.501943,0.998465}v}}{dx}$ we can further expand the equation to get:

$$ \large \require{color} \frac{\gamma -1}{\gamma}\left(-{\color[rgb]{0.814433,0.253157,0.091125}\rho} {\color[rgb]{0.059472,0.501943,0.998465}v}\frac{d{\color[rgb]{0.059472,0.501943,0.998465}v}}{dx} + {\color[rgb]{0.064095,0.501831,0.501977}h}\frac{d{\color[rgb]{0.814433,0.253157,0.091125}\rho}}{dx}\right) + {\color[rgb]{0.814433,0.253157,0.091125}\rho} {\color[rgb]{0.059472,0.501943,0.998465}v} \frac{d{\color[rgb]{0.059472,0.501943,0.998465}v}}{dx} $$

Then we can use that $\require{color}\frac{{\color[rgb]{0.315209,0.728565,0.037706}p}}{{\color[rgb]{0.814433,0.253157,0.091125}\rho}}= \frac{\gamma -1}{\gamma} {\color[rgb]{0.064095,0.501831,0.501977}h}$ to further expand to $\require{color}-\frac{\gamma -1}{\gamma}{\color[rgb]{0.814433,0.253157,0.091125}\rho} {\color[rgb]{0.059472,0.501943,0.998465}v}\frac{d{\color[rgb]{0.059472,0.501943,0.998465}v}}{dx} + \frac{{\color[rgb]{0.315209,0.728565,0.037706}p}}{{\color[rgb]{0.814433,0.253157,0.091125}\rho}}\frac{d{\color[rgb]{0.814433,0.253157,0.091125}\rho}}{dx} + {\color[rgb]{0.814433,0.253157,0.091125}\rho} {\color[rgb]{0.059472,0.501943,0.998465}v} \frac{d{\color[rgb]{0.059472,0.501943,0.998465}v}}{dx}$, simplifying this we get:

$$ \large \require{color}\frac{1-\gamma}{\gamma}{\color[rgb]{0.814433,0.253157,0.091125}\rho} {\color[rgb]{0.059472,0.501943,0.998465}v}\frac{d{\color[rgb]{0.059472,0.501943,0.998465}v}}{dx} + \frac{{\color[rgb]{0.315209,0.728565,0.037706}p}}{{\color[rgb]{0.814433,0.253157,0.091125}\rho}}\frac{d{\color[rgb]{0.814433,0.253157,0.091125}\rho}}{dx} + \frac{\gamma}{\gamma}{\color[rgb]{0.814433,0.253157,0.091125}\rho} {\color[rgb]{0.059472,0.501943,0.998465}v} \frac{d{\color[rgb]{0.059472,0.501943,0.998465}v}}{dx} = \frac{1}{\gamma}{\color[rgb]{0.814433,0.253157,0.091125}\rho} {\color[rgb]{0.059472,0.501943,0.998465}v}\frac{d{\color[rgb]{0.059472,0.501943,0.998465}v}}{dx} + \frac{{\color[rgb]{0.315209,0.728565,0.037706}p}}{{\color[rgb]{0.814433,0.253157,0.091125}\rho}}\frac{d{\color[rgb]{0.814433,0.253157,0.091125}\rho}}{dx} $$

Then using the fact that $\require{color}\frac{d{\color[rgb]{0.814433,0.253157,0.091125}\rho} {\color[rgb]{0.059472,0.501943,0.998465}v}}{dx}=0$ we can get that $\require{color}\frac{d{\color[rgb]{0.814433,0.253157,0.091125}\rho} {\color[rgb]{0.059472,0.501943,0.998465}v}}{dx}={\color[rgb]{0.059472,0.501943,0.998465}v}\frac{d{\color[rgb]{0.814433,0.253157,0.091125}\rho}}{dx} + {\color[rgb]{0.814433,0.253157,0.091125}\rho}\frac{d{\color[rgb]{0.059472,0.501943,0.998465}v}}{dx}=0 \longrightarrow \frac{1}{{\color[rgb]{0.814433,0.253157,0.091125}\rho}}\frac{d{\color[rgb]{0.814433,0.253157,0.091125}\rho}}{dx}=-\frac{1}{{\color[rgb]{0.059472,0.501943,0.998465}v}}\frac{d{\color[rgb]{0.059472,0.501943,0.998465}v}}{dx}$, subsituting this into the equation we get:

$$ \large \require{color}\frac{d{\color[rgb]{0.814433,0.253157,0.091125}\rho} {\color[rgb]{0.059472,0.501943,0.998465}v}}{dx}={\color[rgb]{0.059472,0.501943,0.998465}v}\frac{d{\color[rgb]{0.814433,0.253157,0.091125}\rho}}{dx} + {\color[rgb]{0.814433,0.253157,0.091125}\rho}\frac{d{\color[rgb]{0.059472,0.501943,0.998465}v}}{dx}=0 \longrightarrow \frac{1}{{\color[rgb]{0.814433,0.253157,0.091125}\rho}}\frac{d{\color[rgb]{0.814433,0.253157,0.091125}\rho}}{dx}=-\frac{1}{{\color[rgb]{0.059472,0.501943,0.998465}v}}\frac{d{\color[rgb]{0.059472,0.501943,0.998465}v}}{dx} $$

Lastly by using that $\require{color}{\color[rgb]{0.814433,0.253157,0.091125}\rho} {\color[rgb]{0.059472,0.501943,0.998465}v}^2 = \frac{1}{\gamma}{\color[rgb]{0.315209,0.728565,0.037706}p}{\color[rgb]{0.041893,0.355669,0.727621}M}^2$ we can get that $\require{color} \frac{1}{\gamma}{\color[rgb]{0.814433,0.253157,0.091125}\rho} {\color[rgb]{0.059472,0.501943,0.998465}v} = \frac{{\color[rgb]{0.315209,0.728565,0.037706}p}{\color[rgb]{0.041893,0.355669,0.727621}M}^2}{{\color[rgb]{0.059472,0.501943,0.998465}v}}$, substituting this in we get:

$$ \large \require{color}\frac{{\color[rgb]{0.315209,0.728565,0.037706}p}{\color[rgb]{0.041893,0.355669,0.727621}M}^2}{{\color[rgb]{0.059472,0.501943,0.998465}v}}\frac{d{\color[rgb]{0.059472,0.501943,0.998465}v}}{dx} - \frac{{\color[rgb]{0.315209,0.728565,0.037706}p}}{{\color[rgb]{0.059472,0.501943,0.998465}v}}\frac{d{\color[rgb]{0.059472,0.501943,0.998465}v}}{dx} = \frac{{\color[rgb]{0.315209,0.728565,0.037706}p}({\color[rgb]{0.041893,0.355669,0.727621}M}^2 - 1)}{{\color[rgb]{0.059472,0.501943,0.998465}v}}\frac{d{\color[rgb]{0.059472,0.501943,0.998465}v}}{dx} $$

Putting this back into the original force equation we have that:

$$ \large \require{color}-\frac{4{\color[rgb]{0.990308,0.800015,0.121270}\tau}_w}{D} = \frac{{\color[rgb]{0.315209,0.728565,0.037706}p}({\color[rgb]{0.041893,0.355669,0.727621}M}^2 - 1)}{{\color[rgb]{0.059472,0.501943,0.998465}v}}\frac{d{\color[rgb]{0.059472,0.501943,0.998465}v}}{dx} $$

From dimensional analysis we know that $\require{color}c_f = \frac{{\color[rgb]{0.990308,0.800015,0.121270}\tau}_w}{\frac{1}{2}{\color[rgb]{0.814433,0.253157,0.091125}\rho} {\color[rgb]{0.059472,0.501943,0.998465}v}^2}$, substituting this into the equation gets:

$$ \large \require{color}\frac{4c_f}{D}dx = -\frac{2{\color[rgb]{0.315209,0.728565,0.037706}p}}{{\color[rgb]{0.814433,0.253157,0.091125}\rho} {\color[rgb]{0.059472,0.501943,0.998465}v}^2}({\color[rgb]{0.041893,0.355669,0.727621}M}^2-1)\frac{d{\color[rgb]{0.059472,0.501943,0.998465}v}}{{\color[rgb]{0.059472,0.501943,0.998465}v}} $$

To solve we need to integrate along x from 0 to the final length L and from the initial to final velocities:

$$ \large \require{color}\int_0^L \frac{4c_f}{D}dx = \int_{{\color[rgb]{0.059472,0.501943,0.998465}v}_1}^{{\color[rgb]{0.059472,0.501943,0.998465}v}_2} -\frac{2{\color[rgb]{0.315209,0.728565,0.037706}p}}{{\color[rgb]{0.814433,0.253157,0.091125}\rho} {\color[rgb]{0.059472,0.501943,0.998465}v}^2}({\color[rgb]{0.041893,0.355669,0.727621}M}^2-1)\frac{d{\color[rgb]{0.059472,0.501943,0.998465}v}}{{\color[rgb]{0.059472,0.501943,0.998465}v}} $$

By using that $\require{color}{\color[rgb]{0.814433,0.253157,0.091125}\rho} {\color[rgb]{0.059472,0.501943,0.998465}v}^2 = \frac{}{\gamma}{\color[rgb]{0.315209,0.728565,0.037706}p}{\color[rgb]{0.041893,0.355669,0.727621}M}^2$ we can get that $\require{color}-\frac{2{\color[rgb]{0.315209,0.728565,0.037706}p}}{{\color[rgb]{0.814433,0.253157,0.091125}\rho} {\color[rgb]{0.059472,0.501943,0.998465}v}^2} = -\frac{2}{\gamma {\color[rgb]{0.041893,0.355669,0.727621}M}^2}$, substituting this in we get:

$$ \large \require{color}\int_0^L \frac{4c_f}{D}dx = \int_{{\color[rgb]{0.059472,0.501943,0.998465}v}_1}^{{\color[rgb]{0.059472,0.501943,0.998465}v}_2} -\frac{2}{\gamma {\color[rgb]{0.041893,0.355669,0.727621}M}^2}({\color[rgb]{0.041893,0.355669,0.727621}M}^2-1)\frac{d{\color[rgb]{0.059472,0.501943,0.998465}v}}{{\color[rgb]{0.059472,0.501943,0.998465}v}} $$

We want to only integrate with one variable, in this case ${\color[rgb]{0.041893,0.355669,0.727621}M}^2$ as we want the relationship between Mach and coefficient of friction, and ${\color[rgb]{0.041893,0.355669,0.727621}M}^2$ is easier to integrate over. To do this we start with the equation for stagnation enthalpy $\require{color}{\color[rgb]{0.064095,0.501831,0.501977}h}_0 = {\color[rgb]{0.064095,0.501831,0.501977}h}+\frac{1}{2}{\color[rgb]{0.059472,0.501943,0.998465}v}^2$ and use that $\require{color}{\color[rgb]{0.064095,0.501831,0.501977}h}={\color[rgb]{0.501958,0.501942,0.014744}c_p}{\color[rgb]{0.121820,0.954406,0.966585}T}$, $\require{color}{\color[rgb]{0.059472,0.501943,0.998465}v}={\color[rgb]{0.041893,0.355669,0.727621}M}{\color[rgb]{0.989013,0.435749,0.811750}a}$, and $\require{color}{\color[rgb]{0.989013,0.435749,0.811750}a}=\sqrt{\gamma {\color[rgb]{0.986252,0.007236,0.027423}R}{\color[rgb]{0.121820,0.954406,0.966585}T}}$ to get:

$$ \large \require{color}{\color[rgb]{0.064095,0.501831,0.501977}h}_0 = {\color[rgb]{0.501958,0.501942,0.014744}c_p}{\color[rgb]{0.121820,0.954406,0.966585}T} + \frac{1}{2}{\color[rgb]{0.041893,0.355669,0.727621}M}^2 \gamma {\color[rgb]{0.986252,0.007236,0.027423}R}{\color[rgb]{0.121820,0.954406,0.966585}T} $$

Then using that $\require{color}{\color[rgb]{0.986252,0.007236,0.027423}R}={\color[rgb]{0.501958,0.501942,0.014744}c_p}-c_v$ we can simplify to get:

$$ \large \require{color}{\color[rgb]{0.064095,0.501831,0.501977}h}_0 = {\color[rgb]{0.501958,0.501942,0.014744}c_p}{\color[rgb]{0.121820,0.954406,0.966585}T} + {\color[rgb]{0.501958,0.501942,0.014744}c_p}{\color[rgb]{0.121820,0.954406,0.966585}T}\frac{\gamma -1}{2}{\color[rgb]{0.041893,0.355669,0.727621}M}^2 = {\color[rgb]{0.501958,0.501942,0.014744}c_p}{\color[rgb]{0.121820,0.954406,0.966585}T}(1+\frac{\gamma-1}{2}{\color[rgb]{0.041893,0.355669,0.727621}M}^2) $$

Next we solve for temperature from $\require{color}{\color[rgb]{0.989013,0.435749,0.811750}a}=\sqrt{\gamma {\color[rgb]{0.986252,0.007236,0.027423}R}{\color[rgb]{0.121820,0.954406,0.966585}T}} \longrightarrow {\color[rgb]{0.121820,0.954406,0.966585}T}=\frac{{\color[rgb]{0.989013,0.435749,0.811750}a}^2}{{\color[rgb]{0.986252,0.007236,0.027423}R}\gamma}=\frac{{\color[rgb]{0.059472,0.501943,0.998465}v}^2}{{\color[rgb]{0.041893,0.355669,0.727621}M}^2\gamma {\color[rgb]{0.986252,0.007236,0.027423}R}}$, plugging this in we have:

$$ \large \require{color}{\color[rgb]{0.064095,0.501831,0.501977}h}_0 = {\color[rgb]{0.501958,0.501942,0.014744}c_p}\frac{{\color[rgb]{0.059472,0.501943,0.998465}v}^2}{{\color[rgb]{0.041893,0.355669,0.727621}M}^2\gamma {\color[rgb]{0.986252,0.007236,0.027423}R}}(1+\frac{\gamma-1}{2}{\color[rgb]{0.041893,0.355669,0.727621}M}^2) $$

We need to integrate in order to get a relationship between $d{\color[rgb]{0.059472,0.501943,0.998465}v}$ and $d{\color[rgb]{0.041893,0.355669,0.727621}M}^2$, to make this easier we take the logarithmic of both sides and use logarithmic rules in order to separate the right hand side:

$$ \large \require{color}ln({\color[rgb]{0.064095,0.501831,0.501977}h}_0) = ln\left({\color[rgb]{0.501958,0.501942,0.014744}c_p}\frac{{\color[rgb]{0.059472,0.501943,0.998465}v}^2}{{\color[rgb]{0.041893,0.355669,0.727621}M}^2\gamma {\color[rgb]{0.986252,0.007236,0.027423}R}}\left(1+\frac{\gamma-1}{2}{\color[rgb]{0.041893,0.355669,0.727621}M}^2\right)\right) = \ln\left(\frac{{\color[rgb]{0.501958,0.501942,0.014744}c_p}}{{\color[rgb]{0.986252,0.007236,0.027423}R}\gamma}\right) + ln\left({\color[rgb]{0.059472,0.501943,0.998465}v}^2\right) - ln\left({\color[rgb]{0.041893,0.355669,0.727621}M}^2\right) + ln\left(1+\frac{\gamma -1}{2}{\color[rgb]{0.041893,0.355669,0.727621}M}^2\right) $$

Then we can take the derivative of both sides:

$$ \large \require{color}0 = 0 + \frac{1}{{\color[rgb]{0.059472,0.501943,0.998465}v}^2}(2{\color[rgb]{0.059472,0.501943,0.998465}v}d{\color[rgb]{0.059472,0.501943,0.998465}v}) - \frac{1}{{\color[rgb]{0.041893,0.355669,0.727621}M}^2}d{\color[rgb]{0.041893,0.355669,0.727621}M}^2 + \frac{1}{1+\frac{\gamma-1}{2}{\color[rgb]{0.041893,0.355669,0.727621}M}^2} \frac{\gamma-1}{2}d{\color[rgb]{0.041893,0.355669,0.727621}M}^2 = \frac{d{\color[rgb]{0.059472,0.501943,0.998465}v}}{{\color[rgb]{0.059472,0.501943,0.998465}v}} - \frac{d{\color[rgb]{0.041893,0.355669,0.727621}M}^2}{{\color[rgb]{0.041893,0.355669,0.727621}M}^2\left(1+\frac{\gamma -1}{2}{\color[rgb]{0.041893,0.355669,0.727621}M}^2\right)} $$

Solving for $\frac{d{\color[rgb]{0.059472,0.501943,0.998465}v}}{{\color[rgb]{0.059472,0.501943,0.998465}v}}$ we get that:

$$ \large \require{color}\frac{d{\color[rgb]{0.059472,0.501943,0.998465}v}}{{\color[rgb]{0.059472,0.501943,0.998465}v}} = \frac{d{\color[rgb]{0.041893,0.355669,0.727621}M}^2}{2{\color[rgb]{0.041893,0.355669,0.727621}M}^2\left(1+\frac{\gamma-1}{2}{\color[rgb]{0.041893,0.355669,0.727621}M}^2\right)} $$

Then by substituting this into the integral we have:

$$ \large \int_0^L \frac{4c_f}{D}dx = -\int_{{\color[rgb]{0.041893,0.355669,0.727621}M}_1}^{{\color[rgb]{0.041893,0.355669,0.727621}M}_2} \frac{1}{\gamma {\color[rgb]{0.041893,0.355669,0.727621}M}^2} (1-{\color[rgb]{0.041893,0.355669,0.727621}M}^2) \frac{d{\color[rgb]{0.041893,0.355669,0.727621}M}^2}{{\color[rgb]{0.041893,0.355669,0.727621}M}^2\left(1+\frac{(\gamma-1)}{2}{\color[rgb]{0.041893,0.355669,0.727621}M}^2\right)} $$

Integrating with ${\color[rgb]{0.041893,0.355669,0.727621}M}_2 =1$ we get that:

$$ \large \require{color}\frac{4c_{f}}{D}L = \frac{1 - {\color[rgb]{0.041893,0.355669,0.727621}M}_{1}^{2}}{\gamma {\color[rgb]{0.041893,0.355669,0.727621}M}_1^2} + \frac{\gamma +1}{2\gamma} \ ln\left(\frac{\frac{\gamma +1}{2}{\color[rgb]{0.041893,0.355669,0.727621}M}_1^2}{1+\frac{\gamma -1}{2}{\color[rgb]{0.041893,0.355669,0.727621}M}_1^2}\right) $$